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ABSTRACT 

We present a comprehensive set of axisymmetric, time-dependent simulations of 
jets from Keplerian disks whose mass loading as a function of disk radius is systemat- 
ically changed. For a reasonable model for the density structure and injection speed 
of the underlying accretion disk, mass loading is determined by the radial structure 
of the disk's magnetic field structure. We vary this structure by using four different 
magnetic field configurations, ranging from the "potential" configuration (Ouyed & 
Pudritz 1997), to the increasingly more steeply falling Blandford & Payne (1982) and 
Pelletier & Pudritz (1992) models, and ending with a quite steeply raked configuration 
that bears similarities to the Shu X-wind model. We find that the radial distribution 
of the mass load has a profound effect on both the rotational profile of the underlying 
jet as well as the degree of collimation of its outflow velocity and magnetic field lines. 
These four models have systematic differences in the power-law rotation profiles of jet 
material far from the source; i>^( r ) oc r a ranging over —0.46 ^ a ^ —0.76. We show 
analytically, and confirm by our simulations, that the collimation of a jet depends on 
its radial current distribution, which in turn is prescribed by the mass load. Models 
with steeply descending mass loads have strong toroidal fields, and these collimate 
to cylinders (this includes the Ouyed-Pudritz and Blandford-Payne outflows). On the 
other hand, the more gradually descending mass load profiles (the PP92 and monopo- 
lar distributions) have weaker toroidal fields, and these result in wide-angle outflows 
with parabolic collimation. We also present detailed structural information about jets 
such as their radial profiles of jet density, toroidal magnetic field, and poloidal jet 
speed, as well as an analysis of the bulk energetics of our different simulations. Our 
results are in excellent agreement with the predictions of asymptotic collimation for 
axisymmetric, stationary jets. 

Key words: accretion: accretion disks - stars: formation - stars: pre-main-sequence 
- ISM: jets and outflows - MHD 



1 INTRODUCTION 

Hydromagnetic disk winds provide what is arguably the 
most comprehensive quantitative picture for the origin of 
astrophysical jets in protostellar systems. An extensive 
body of theoretical work as well as a wide range of nu- 
merical simulations, using very different codes, robustly 
demonstrate that a centrifugally driven wind from an ac- 
cretion disk can efficiently extract its angular momen- 
tum and gravitational potential energy. Disk winds pro- 
vide a potentially universal mechanism for the origins of 
jets that are observed in both lo w and high mass protostar s 
(eg. reviews: iRicher et all (|2000D : iReipurth fc Ballvl feOOll) : 
IZhand lEislhlShepherdf (2005)) as well many other systems 



such as micr oquasars (eg . iMirabel fc Rodriguez! <| l998f0 and 
quasars (e.g. lBeeelman. Blan^ford^^R^eain^84|) ) because 
jet properties - such as their terminal velocities and wind 
mass-loss rates - naturally scale with the depth of the gravi- 
tational potential well of t he central obje c t and the accretion 
rate i n to it (eg . revie ws; ISpruiti jl99rJ); IKonigl fc Pudrit j 
j2000D : lPudrit3 i2003l) ; fPudritz fc Baneried J2005D 1. 

The theory has recently received considerable empiri- 
cal support by HST observations of the rotati on of several 
jets a s sociated with T-Tauri stars (TTS) (eg. ICoffev et alJ 
j2004 : lBacciotti et alJ j2000ll2002D ). High spectral and spa- 
tial resolution observations have directly measured the ro- 
tation of these jets. The results suggest that the observed 
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jets might have arised from a region of roughly 1 AU in 
size in the disk. Moreover, the angular momentum that jets 
are inferred to carry from this measurement, is a significant 
fraction that must be carried off to drive the observed ac- 
cretion rates in the underlying disks. These results seem to 
agre e with the theoretical predictions of disk wind theory 
(eg.. lAnderson et all (120031) ). 

The observations of jets provoke several important ques- 
tions about the role of disk winds. How do the underlying 
accretion disks control jet collimation, even to very large 
scales? And what determines how rapidly the observed jets 
spin and the amount of angular momentum that they can 
extract from the disk? We show in this paper that both of 
these important properties of jets are aspects of the toroidal 
dynamics of the jet and are controlled by the same physical 
process, namely, by the mass loading of the disk wind field 
lines that occurs near the surface of the disk. 

Consider first the issue of jet collimation. The basic 
mechanism has been understood, at least in principle, for 
some time, ft involves the action of a hoop stress that arises 
from the generation of a toroidal magnetic field in the ro- 
tating, magnetized flow, which pinches the flow towards the 
outflow axis. The detailed analysis for the collimation of out- 
flows for steady-state jets is extremely difficult since it in- 
volves the solution of the highly nonlinear, Grad-Shafranov 
equation. Progress has been limited to studies of the asymp- 
totic properties of steady state jets l|Hevvaerts fc Normanl 
1989) (hen cefor th HN 89), o r s pecial self-similar models (eg. 
bstrikerl Jl99St): O dl995l): iTrussoni. Tsinganos fc Sautvl 
Jl997l) : ICasse fc Ferreiral feOOOl) ). The general collimation 
properties of time-dependent jets have, so far, defied gen- 
eral theoretical analysis. 

Nature probably produces jets with a range of intrinsic 
collimations, as suggested by the observations of molecular 
outflows. Models for observed outflows fall into two general 
categories: the jet-driven bow shock picture, and a wind- 
driven shell picture in which the molecular gas is driven 
by an underlying wide-angle wind component (eg. review 
ICabrit. Raga. fc Guethl lll997l) ). A survey of molecular out- 
flows bv lLee et alJ (I2000T) found that both mechanisms are 
needed in order to explain the full set of systems observed. 

Can the mechanism of jet collimation readily ac- 
commodate both possibilities? The simulations of disk 
winds explore a range of initial magnetic configurations 
that either origina te on the surface of a star and c onnec t 
with the disk (eg. iHjarashj^Jjhibj^^ 1199^); 
iGoodson Winglee. fc Bohml dl997l) : iFendt fc Elstnerl 
J1999D: iKeppensOOl fcOOOD) with in the disk (eg., 
Ivon Rekowski fc Brande nburg| (|2004j)), or that threa d 
only the d isk it sel f (eg. lOuyed, Pudritz. fc Stond Jl997f 



Tomisakal <ll998l): [Kudoh. Matsumoto. fc Shibatal 120021) : 
Fendt fc Cemeliid i2002f) ). In the former case, twisting of 



the magnetic field that initially threads the disk beyond the 
co-rotation radius inflates and then disconnects the field 
from the star to produce a steeply decreasing, disk-field. 
In the latter case of purely disk-threaded field, simulations 
have a variety of initial magnetic geom etries: from the 
split monopole of iRomanova et alJ fll997TI . similar to the 



magnetic geometry of Blandford & Payne (1982; BP82) is 
intermediate between these two extreme regimes with an 
initial poloidal magnetic field on the surface of the disk 
that falls off with disk radius as Bbps,2(t ) oc ro '^ A - 

All of these simulations show that outflows are centrifu- 
gally driven but that the degree to which they collimate 
varies. Although the steeply raked, monopolar-like, mag- 
netic field distributions fail to collimate to the outflow axis, 
they are still candidates for jet models. This is because the 
densest part of the outflow is directed along those more po- 
lar lines giving the appearance of "collimated jets" (even 
wind s from purely dipole stellar field pr o duce this struc- 
ture:ISakuraaJl985l) : lMatt fc Balicld(l2004h : lMatt fc Pudritj 



highly concentrated mag netic field of the IShu et al. ( 2000 
X-wind configuration, to the far more gradual declining or 
constant magnetic fields featured in the potential solutions 
investigated by OPS97. The familiar, self-similar disk wind 
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We show in this paper that it is the mass loading of 
a jet - defined as the mass outflow rate per unit magnetic 
flux - that determines how well jets are collimated. In earlier 
theory and numerical papers we found that low mass loads 
for jets lead to rapid, episodic behaviour while more heav- 
ily mass-loaded systems tend to achieve stationary outflow 
configurations (OPS97, Ouyed & Pudritz 1997b (OP97b), 
Ouyed & Pudritz 1999 (OP99)). This has al so been exam- 
ined in recent work bv lAnderson et al.1 (120051) . 

Regarding the second major issue - angular momentum 
transport in the jet - we show that this is related to an- 
gular momentum extraction by the disk wind torque upon 
the disk. In steady state, the wind torque depends on the 
strength of the toroidal magnetic field near the disk surface 
and this too, we show, depends on how the wind is mass 
loaded at the disk surface. 

We present a series of numerical simulations of jets us- 
ing a much broader range of initial magnetic configurations 
than employed in our earlier papers. We vary the mass per 
unit magnetic flux that threads the disk (the mass load- 
ing function) by employing a wide range of initial magnetic 
field distributions across our disk model. Our initial mag- 
netic configurations range from rather well collimated ini- 
tial fields (such as the "potential configuration" of OP97a) 
and the BP82 configuration, to the initially less well colli- 
mated configuration of Pelletier-Pudritz (1992; PP92). The 
final configuration in our simulation suite features an even 
more steeply declining threading field with disk radius, and 
is even more open. 

One of our major findings is that there are indeed two 
different kinds of collimation that can be achieved by disk 
winds with different mass loading profiles: cylindrical colli- 
mation and a wide-angle outflow. We provide detailed nu- 
merical data that show the radial structure and collimation 
of a variety of physical variables, including all components 
of the velocity and magnetic field of the jet, as well as its 
density profile. We show that this agrees with the predic- 
tions that jet collimation is related to the radial structure 
of the current in these jets, as was first demonstrated ana- 
lytically by Heyvaerts & Norman (1989, HN89). In fact, the 
wide-angle flows we find correspond quite well to the case of 
parabolic collimation predicted in this theory. 

We also show that the mass loading of jets affects their 
radial angular rotation profiles. In particular, we find that 
the radial distribution of jet rotation is robust and follows 
radial power-law relations, oc r a , whose index a varies 
little from model to model: typically in the range —0.46 ^ 
a ^ —0.76 depending on the mass-load profile. These radial 
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rotation profiles can in principle, be used to discriminate 
between different disk wind models because the distribution 
of angular momentum in the wind reflects on how it is mass 
loaded by the accretion disk at its base. 

We begin this paper by first briefly summarizing some 
recent important observational results on the relation be- 
tween accretion disks and jets from YSOs (§2) and follow- 
ing this up with a brief review of the basic aspects of hydro- 
magnetic disk wind theory needed to interpret our numerical 
results (§3). We then outline our generalized disk wind mod- 
els which include simulations of BP82 and PP92, and other 
mass loading profiles, in §4. We analyze our results in §5, 
and conclude in §6. 



2 OBSERVATIONS OF JET STRUCTURE AND 
DYNAMICS 

Of the many important contributions on the issue of jet 
structure and dynamics, we focus on just a few recent ad- 
vances that theories and simulations must address. As an 
example, recent detailed high resolution imaging and spec- 
troscopic studies of jets have revealed their internal dynami- 
cal structure in sufficient detail to provide useful constraints 
on jet theories and numerical simulations. High resolution 
images of H a and [01] emission in the DG Tau jet taken in 
various velocity ranges from 50 - 450 km s" 1 in 125 km s" 1 
wide velocity bins feacciotti et al.ll2002l) reveal an "onion- 
like" velocity structure wherein there is a continuous vari- 
ation in the velocity of the jet, the very highest being the 
densest and nearest the jet axis, and the lowest farther away. 
The lower speed flow is farther from the axis and appears 
to broaden out and become less collimated with increasing 
radius in the jet. The jet can be traced to within 15 AU 
(0.1") of the central source. While earlier papers suggested 
that jets may have two distinct velocity components (fast 
and hig hly collimated, versus slow er and wide- angled com- 
ponents lkwan fc Tademarul |l988)), the recent data suggest 
a continuously varying velocity structure. 

A second important observed property of jets is that 
they rotate. For DG Tau as an example, the flow appears to 
be rota ting clockwise when vie wed from the flow towards the 
source iBacciotti et alJ|2003allbT) . The rotation speed is 6 - 
15 km s , depending on the position in the flow within 110 
AU from the source. These authors find agreement between 
these measured rotation speeds with those expected in a 
centrifugally driven hydromagnetic wind. Using the scalings 
for disk wind speeds and angular momentum transport rates 
outlined in §3.1 and §3.2, this outflow is estimated to be 
carrying up to 60 % of the angular momentum loss from the 
disk at its footpoint that is needed to drive the observed 
accretion rate. 

Finally, as has been noted many times in the literature, 
there is a strong correlation between the observed mass ac- 
cretion rate through protostellar disks and the mass loss 
rates carried by their jets. Optical jets associated with TTS 
have typical mass loss rates that are consistently about a 
tenth of the underly i ng dis k accretion rate. As an example, 
lHartmann fc Calved dl995l) found that for a sample of 31 
TTS, jet mass loss rates calculated from the forbidden lines 
in the range 3 x 10~ 7 < Mj < 10" 10 M o yr _1 with a median 
value of 1O" 9 M yr _1 (eg. review Calvet 2003) . These mass 



loss rates are a tenth of the underlying disk accretion rate 
in the system; 

Mj/M a ~ 0.1. (1) 

Theory predicts that this is a reflection of the efficient angu- 
lar momentum extraction from accretion disks by the wind 
torque. 



3 DISK WIND THEORY: JET STRUCTURE 
AND THE ROLE OF MASS LOADING 

The theory of hydromagnetic winds originated with t he early 
work o n winds from ro t ating magnetized stars (eg. iMestell 
Jl968ft : IWeber fc Davis! <196m These papers developed 1- 
D, axisymmetric models of hydromagnetic flows from ro- 
tating stars. The application of this idea to self-similar ac- 
cretion disks was first carried out in the seminal paper by 
BP82. 

Consider the simplest possible description of magne- 
tized, rotating, conducting gas of density p that is threaded 
by a large-scale field. There are 4 equations of stationary, 
ideal MHD that describe such a system, namely: the conser- 
vation of mass (continuity equation) ; the equation of motion 
for gas undergoing a pressure gradient force (from the pres- 
sure p), as well as a Lorentz force (from the field B) - all 
within in the gravitational field of a central object (whose 
gravitational potential is <j>); the induction equation for the 
evolution of the magnetic field in the moving gas; and fi- 
nally, the conservation of magnetic flux. These equations 
were written down by Chandrasekhar, Mestel, and many 
others: 

V.(pv) = (2) 



pv.Vv = -Vp-pV(j>+ — (V x B) x B (3) 

47T 

V x (v x B) = (4) 

V.B = (5) 

It is convenient to decompose all magnetic fields and velocity 
fields into poloidal and toroidal field components B = B p + 
B^e^ and v = v p + v^e^. 

The most basic conservation laws that these equations 
express are the conservation of mass and magnetic flux 
(equations 2 and 5). By comparing the form of these two 
equations, we see that the poloidal field and the poloidal 
mass flux along the field line must be proportional to one 
another along a field line. Defining a function k that is a 
constant along a magnetic field line and which we will call 
the "mass load" of the wind, then; 

pv p = fcBp (6) 

This mass load function can be cast in a much more 
revealing way by noting that the wind mass loss rate passing 
through an annular section of the flow of area dA through 
the flow is dM w = pv p dA, while the amount of poloidal 
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magnetic flux through this same annulus is d$ = B p dA. 
Thus, the mass load per unit time, per unit magnetic flux of 
the wind that is preserved along each streamline along the 
flow emanating from the rotor (a disk in this case) is 



B„ 



dM w 



(7) 



The ultimate source for the mass load is, of course, the 
accretion disk. While there may be several complicated pro- 
cesses in the disk that prescribe the function k, the station- 
ary theory summarized in equation (7) states that this mass 
load all comes down to the playoff between the gas density, 
the injection speed, and the magnetic field strength across 
the surface of the accretion disk. 

The strength of the toroidal magnetic field in the jet 
plays an important role in jet collimation. We calculate the 
strength of this field component at any point in a stationary, 
axisymmetric jet by using the poloidal component of the 
induction equation which can be written as, 



o ) 



(8) 



This result shows that the strength of the toroidal field at 
any point in the flow depends on three things: the relative 
shear of the flow between the point in question and the foot- 
point of the field line; the density of the gas (the denser it 
is the greater the inertia and hence the greater the toroidal 
field); and finally the mass load. The larger the value of k 
along a field line, the smaller is the value of the toroidal 
field, and hence the more diminished role that it plays in 
collimating the outflow. All else being equal, we expect to 
have stronger toroidal fields in the portions of an outflow 
that are fed with smaller mass loads. From this it is clear 
that the radial distribution of the mass load k in the jet will 
play a significant role in controlling its collimation. 

Now let us turn our attention to the radial rotational 
profile of a hydromagnetic jet. We derive this from the con- 
servation of angular momentum equation for axisymmetric, 
stationary flows. This is described by the <j> component of 
equation (3), which yields; 



pv p .V(rt>0) 



f*.V(rB,) 



(9) 



Let us apply this equation to two different situations - (i) 
the transport of angular momentum along a field line and 
(ii) the transport of angular momentum out of a disk, by an 
outflow (which is analyzed in the following subsection). 

By applying equation (6) to equation (9) the angular 
momentum equation reduces to 

.V ( ™,-^)=0 (10) 



B, : 



This equation says that there is a conserved quantity, 

rB^ 



I — rv^ 



4nk' 



(11) 



along any given field-line in the flow. This is the total an- 
gular momentum per unit mass. The form for I reveals that 
the total angular momentum is carried by both the rotating 
gas (first term) as well by the twisted field (second term). 
Once again, we see that the mass load plays a significant 
role - here by controlling the relative amount of angular mo- 
mentum that is carried by the rotation of material in the jet 



as compared to that carried by the jet's twisted magnetic 
field. 

The conservation of angular momentum given above, 
can be rearranged to find the rotational profile of the jet; 



Im — r 0, o 
m 2 — 1 



(12) 



where the so-called Alfven Mach number of the flow is m — 
Vp/vA and va = B p / \f\A'Kp) is the Alfven speed of the flow. 
The Alfven surface is that point r — ta at each field line in 
the wind where m — 1. The flow may be regarded as kept 
in co-rotation with the disk until this point is reached. 

The actual value that I takes along any field line is de- 
termined by the Alfven surface where m = 1; 



1(a) 



Q. r\ 



(13) 



For a field line starting at a point r on the rotor (disk in 
our case), the Alfven radius is rA(r ) and constitutes a lever 
arm for the flow. 

Finally, we note that the conservation of energy in a 
stationary jet yields a simple relation between the terminal 
speed of a jet, and the Alfven radius: 

Voo ~ 2 1/2 (r A /r )vK{r ) (14) 

where VK(r ) is the Kepler speed at the base of the field line 
that threads the disk at radius r . By combi ning the conse r- 
vation laws, one may show that (eg. PP92. ISpruitl <199eft 1: 



ta/to = [(B 2 /2Trpv p Q r ] 1/:i = ^ 



1/3 



(15) 



The dimensionless parameter n - which is not a constant 
along a given field line - ha s been used as di f ferent measure 
of the mass load of a jet bv I Anderson et ail i2005h . 

In general, the conservation laws show that the mass 
load plays a fundamental role in determining both the col- 
limation as well as the rotational profile of the jet. Even 
though direct communication of the jet with the accretion 
disk is cut off beyond the fast magnetosonic surface of the 
outflow (which is rather close to the disk) - nevertheless the 
mass load function is transported from the accretion disk to 
every point in the outflow. This is how the accretion disk can 
control the jet's dynamics even far from the source. 



3.1 Jet rotation: angular momentum extraction 
from the disk 

Let us now carry out our second application of the angular 
momentum equation (9), this time to find the torque that 
is exerted upon a thin accretion disk by an external field B. 
Any vertical flow in the thin disk is negligible so only the 
v r contribution matters, and the rotation speed v$ is the 
Kepler speed if disks are thin. After vertically integrating 
the equation, noting that the disk accretion rate is M a = 
— 2nTiV r r . The result is the angular momentum equation for 
the accretion disk under the action of an external magnetic 
torque; 



, v d{r v ) _ 2 

Ma j = — r o D$tiz\ ro ,H 

ar o 



(16) 



This result shows that angular momentum is extracted out 
of disks threaded by magnetic fields. By solving for rB$ = 
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k(rv$ — I) and relations (7) and (11) for k and I the disk 
angular momentum equation (14) can be cast into its most 
fundamental form; 



d{n o r 2 ) _ dM u 



dr a 



dr a 



n„r\.(l - {ro/r A ) 2 ) (17) 



As is now well understood, this fundamental equation re- 
veals the crucial link between the mass outflow in the wind, 
and the mass accretion rate through the disk; 



M a = {r A /r fMu. 



(18) 



The magnetic field forces gas to co-rotate with the un- 
derlying disk out to a distance of the order ta which is 
why these winds can be so efficient at extracting the angu- 
lar momentum of the disk. Later in the paper we will show 
that the typical Alfven lever arm values in the simulations 
is ta/to — 2 — 3 for a broad range of models. The mean 
value of the lever arm allows one to calculate the ratio of 
the wind mass outflow rate compared to the disk accretion 
rate, M w /M a ~ 0.1 - 0.3. 



3.2 Two possible solutions for jet collimation: 

cylindrically collimated vs "wide-angle" flows 

In the standard picture of hydromagnetic winds, collima- 
tion of an outflow occurs because of the increasing toroidal 
magnetic field in the flow as one moves through its various 
critical points. From equation (8) one deduces that at the 
Alfven surface B^, ~ B p , while in the far field (assuming that 
rj >> ta) this ratio is of the order B^/B p ~ r/rA- There- 
fore, collimation can in principle be achieved by the tension 
force associated with the toroidal field which leads to a ra- 
dially inwards directed component of the Lorentz force (or 
"z-pinch"); F Lore „t^ r ~ J^B^. 

The radial structure of the outflow is found from the 
condition of force balance perpendicular to the field lines 
in the flow, which is known as the Grad-Shafranov equa- 
tion. This is a very complicated, non-lin ear equation and n o 
general solutions are known (see review. iHevvaerta (|2003)). 
Because of the mathematical difficulties, analytic studies 
have been dominated by simplified approaches which focus 
mainly on f inding stationary, self-s i milar soluti ons of the flow 
(eg. BP82. pTsinganos fc Trussonil dl990lk [LI lll995D ; iKonigll 
jl989|hlOstriked il998MCasse fc Ferreiral i2000h '). or on ana- 
lyzing the asymptotic limits of the general equations (HN89, 
PP92). 

In HN89 it was shown that two types of solution are 
possible depending upon the asymptotic behaviour of the 
total current intensity in the jet; 



I = 2tt 



/ Jz(r',z')dr' = (c/2)rB 

J 



(19) 



where J = (c/47r)V x B is the current density. In the limit 
that / — > as r — * oo, the field lines are paraboloids which 
fill space. On the other hand, if the current is finite in this 
limit, then the flow is collimated to cylinders. The collima- 
tion of a jet therefore depends upon its current distribution 
- and hence on the radial distribution of its toroidal field - 
which, as we saw earlier, depends on the mass load. Mass 
loading therefore determines the asymptotic behaviour of 
jets! 



In order to provide a context for the simulation re- 
sults to come, it is useful to briefly overview the self-similar 
model of BP82 in which the disk accretion rate treated 
as a constant. The only velocity in the self-similar prob- 
lem is the Kepler speed of the disk. Therefore, one antic- 
ipates that the various velocities in the problem scale as 
va oc c a oc v r oc Voo oc v k where one has the Alfven, sound, 
radial inflow, terminal wind speed, and Kepler speed all be- 
ing proportional to one another. Now, the hydrostatic bal- 
ance condition for thin disk gives H(r)/r = c 3 /vk- There- 
fore, the self-similar scaling c s oc vk implies that H oc r; the 
disk is wedge-like in its spatial structure. This also implies 
that the disk temperature follows the virial scaling T oc r _1 . 
The density scaling for the disk follows from the definition 
of the disk accretion rate, M a = 2tt(2H p)v r .r, which given 
the scaling of the radial inflow velocity v r oc vk and the disk 
scale height relation implies that p oc r _3//2 . The wind mass 
loss rate in the wind follows since the mass loss per unit area 
along the surface of the disk is dM w jdA — p V( w , ) oc r^ 2 , 
which implies that M m (r ) oc lnr . 

Note that the scaling of the disk Alfven speed va oc vk, 
together with the density implies a particular scaling of the 
disk poloidal field, B p oc r~ 5/4 . In the BP82 model, the 
toroidal field on the disk should therefore scale as B^ /B p = 
const, so that B^ oc r -5 ^ 4 . This presents a bit of a dilemma 
however because the lowest energy toroidal field would be 
expected to be a force-free condition which has a different 
scaling B^ oc r" 1 . The Alfven surface in this model is then 
Voo/vk oc ta/to = const; ie, a cone. Finally, given the scal- 
ing for the disk magnetic field, we can now calculate the 
mass load of the wind in the BP82 model; 

k , B ps2 oc r" 3/4 . (20) 

Jets need not be self-similar structures however. In sta- 
tionary states, it is perhaps more natural to think of them 
as having attained minimum energy configurations instead. 
We have seen above that for the magnetic configuration in 
which the toroidal field scales as B^ oc rj 1 , the current in- 
tensity is finite everywhere. The BP82 self-similar solution 
leads to diverging currents which corresponds to a higher 
energy state for the magnetic configuration. The disk-wind 
solutions elucidated by PP92 were designed to characterize 
non self-similar, minimum energy configurations. 

We briefly remind the reader of the motivation for 
the PP92 solution by first writing down the general Grad- 
Shafranov (GS) equation for stationary, axisymmetric flow, 
for the surfaces of magnetic flux a(r, z) in the flow. Magnetic 
field lines are restricted to surfaces of constant magnetic flux, 
a(r, z) = const, where the the poloidal magnetic field follows 
from the relation B p = (l/r)V(a)xc5. The GS equation can 
be written in a suggestive form (see PP92): 

J-Va) = \(a)ij(a)a (21) 



where the function 77(a) ~ const and the function 



A(a) 



lr\ AnpA (a)ngrj ( dlnr A \ , , 

— ' A = i^J (22) 



is a constant on a surface of constant magnetic flux and 
depends upon the density at the Alfven point, pA- Solving 
this last relation for the Alfven radius, we see that ta oc 
A 1 /* 3 which is such a weak dependence that solutions with 
A = const were sought (PP92). 
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The accretion disk imposes an important set of bound- 
ary conditions on the GS equation that includes the value 
of the magnetic flux function across its surface; a(r, z) = 
a (r , z ). Without loss of generality, one may assume that 
the flux takes on a power-law form on the disk; 



(23) 



so that the poloidal field threading the disk, which is derived 
from this flux, becomes 



B z (r o ,0) = br" a 



(24) 



The general form of the GS equation given above can 
be exploited by using the Ansatz that A = const, as shown 
by PP92. In this regime, we find many of the interesting 
solutions including the self-similar case of BP82. The GS 
equations show that the Alfven surface takes the form: 



r A /r oc a [4( " +1 »- 31/[2 "' +1)1 oc r ^ +(1/2) 



(25) 



while the current takes the form 

I(r,z) ex a [i-2(M+D]/[2(M+i)] K r -M-d/2). (2 6) 

This result for the radial current distribution contains 
the key insight about the collimation of jets from accre- 
tion disks. We note that there are two distinct categories of 
mass loads that by the HN89 theorem, will have different 
collimation properties: For magnetic field distributions with 
— 1/2 > /x; the radial current profile diverges, 



Zim r ^ cx) J(r) 



(27) 



so that wide-angle flows are to be expected. On the other 
hand, distributions with —1/2 < /j, have currents that vanish 
at infinite radius, 



lim r 



./(r)-»0 



(28) 



so that the flow should collimate to cylinders. 

The particular state investigated by PP92 is the unique 
case in which the current is finite, but does not diverge; 
H = —1/2, whose radial current profile is 



Ipp92 = const, 



(29) 



and poloidal field profile in the disk is B p oc r a 3 ^ 2 . The 
Alfven surface for the PP92 model scales as (rA/r )pp92 oc 

— 1/2 

r a , while the terminal speed in this solution has the be- 
haviour Voa oc r" 1 . Finally, the mass loss rate in the wind 
for the PP92 solution is M w oc [(r /ri) — 1], where ri is the 
inner edge of the disk. The mass loading of this model wind 
scales as 



k 0: pp92 oc r a 



1/2 



(30) 



How does an accretion disk establish the radial profile 
of the outflow mass load function, k(r )? The density at 
the disk surface is prescribed by the hydrostatic balance 
condition with the base of the disk corona above it (the 
initial condition we have used in our own simulations, see 
§4). In reality, we expect that the jet originates from some 
region in the disk corona so that the principle of hydrostatic 
balance between the disk and the base of the corona should 
remain reasonably secure. The second factor in the mass 
load is the injection speed of material into the base of the 
outflow, Vi n j . This speed must be subsonic in magnitude. We 
consider that the most natural scaling of this injection speed 



is with the local Keplerian velocity at each radius of the disk; 
typically Vixxj — 10 vk at any point on the disk. These two 
physically reasonable conditions imply that the main factor 
in determining the mass loading of the disk wind is actually 
the radial distribution of the threading magnetic field across 
the disk. Assuming a power-law scaling of this poloidal field 
component given above, we find that the mass load function 
at the disk surface follows a power law relation; 



k(r a ) = p Vp, /Bp, ococ r 1_ 



(31) 



How do the mass loads of the different wind models 
compare? For the BP82 and PP92 mass-load functions as an 
example: the PP92 mass loading function is a less steeply 
raked function of disk radius than the BP82 self-similar case 
(wherein fcsps2 oc r^ 3 ^ 4 ). The associated current profile for 
the BP82 solution (which can also be found from equations 
(19) and (20) above (with A = const) and pertains to the 
case jj, — —1/4. The current profile in the BP82 theory is 
therefore 



Ibps2 oc r 



1/4 



(32) 



The essential point is that whereas the BP82 current 
I{r)pps2 — > as r — » oo, the Pelletier-Pudritz has a finite 
current in this limit I(r)ppg2 7^ 0. This implies that yet 
steeper mass load profiles and more general ly, the monopole- 
like m agnetic configuratio n of models like iRomanova et alJ 
Jl997l) . and the X-wind of IShu et all fcOOOI) . should diverge 
in this limit. The mass load of these latter models produces a 
toroidal magnetic field that is strongly concentrated towards 
the outflow axis, with little left in the wide-reaches of the 
outflow to effectively collimate the outflow to the axis. 

Thus, we have shown that by using the HN89 limit to- 
gether with the theory of mass loading shown above, that the 
PP92 case should provide a dividing line between jets that 
ultimate collimate to cylinders, and jets whose magnetic and 
streamlines are parabolic and therefore "wide-angled". This 
result implies that disk winds can in principle provide a 
wide-angle pressure that can drive molecular outflows, and 
that this is essentially determined by how the accretion disk 
mass-loads the jet. 



4 NUMERICAL SIMULATIONS - INITIAL 
CONDITIONS 

Our numerical approach in this paper is a gener- 
alization of the method used in our earlier papers 
(OPS97, OP97a,b; OP99). Conditions on the underly- 
ing accretion disk remains constant in time and pro- 
vide a set of fixed boundary conditions for the out- 
flow simulatio n. Several grou p s ha ve exploited this 
appro a ch (eg. lUstvueova et all dl998f): IRomanova et al 



Meier et all (Il99j) :_l Krasn opolskv. Li. fc Blandford 



Fendt fc Cemeli 



The published simula- 
tions differ in their assumed initial conditions, such as 
the magnetic field distribution on the disk, the plasma /3 
(= Pgas/fmag) above the disk surfaces, the state of the initial 
disk corona, and the handling of the gravity of the central 
star. 

We used the ZEUS 2-D code which is arguably the 
best documented and utilized MHD code in the literature 
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Figure 1. Left panels: initial magnetic field configurations for winds with initial potential (/i = 0), Blandford & Payne (fi = —0.25), 
Pelletier &; Pudritz (fi = —0.5), and 7Q (fi = —0.75) configurations shown. Right panels: final magnetic field configurations (at t = 400) 
for each case, with Alfven points (filled circles) and FM points (stars) marked. Note the more open magnetic - and stream line structures 
as the initial magnetic configuration (value of fi) changes. 



{Stone fc Normanl ijlQfll I1994K It is an explicit, finite dif- 
ference code that runs on a staggered grid. The evolution of 
the magnetic field is followed by the method of constrained 
transport. In this approach, if V.B = holds for the initial 
magnetic configuration, then it remains so for all later times 
to machine accuracy. The obvious way of securing this con- 
dition is to use an initial vector potential A(r, z, t = 0) that 
describes the desired initial magnetic field at every point in 
the computational domain. 

The accretion disk in all of our work so far is chosen 
to be initially surrounded by a polytropic corona (7 = 5/3) 
that is in hydrostatic balance with the central object. We 
do this because it is the equation of state used in BP82 
which provides an invaluable analytic and numerical solution 



to which we can compare our simulations. The accretion 
disk at the base of the corona is given a density profile that 
it maintains at all times in the simulation, since the disk 
boundary conditions are applied to the "ghost zones" and 
are not part of the computational domain. It is chosen so 
that the disk surface is in pressure balance with the corona 
above it (ie. p oc r^ 3,72 for a 7 = 5/3 model - as in BP82). 
This hydrostatic state has a simple analytic solution which 
was used as the initial state for all of our simulations. 
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(fi= -0.25, z=0.08) 
(a=-0.778035,b=0. 209857) 




0.0005 - 
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(fi= 0.0, z = 72.0801) 
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(fi= -0.25, z = 72.0801) 

(a=-0. 33821 3,b=0. 00383 772) 




(a) Radial dependence of the mass loading parameter, k at the 
disk surface. The first two cases fit the predicted behaviour; k ex 
r -1- ^. (A power law form ip = br a is the fit used for any physical 
quantity tp). 



(b) Radial dependence of density for each configuration at high 
z. The peak density (which defines the jet) in the first two cases 
is about an order of magnitude greater than latter two, and is 
collimated closer to the disk axis. 



(fj.= 0.0, z = 72.0801) 











(a=-0.77918C 


b=4.91S95) 



(li= -0.25, z = 72.0801) 



(li- -0.5, z = 72.0801) 





0.02 - 
0.00 - 



()i= 0.0, z = 72.0801) 
(a=-0.6231B2, b = -0 1541 St j 



(fi- -0.25, z-72.0801) 




(c) Radial dependence of the poloidal velocity at high z. The 
collimated flows show higher velocities closer to the disk axis, 
whereas the non-collimated solutions have almost the opposite 
behaviour, indicative of a wide-angle wind. 



(d) Radial dependence of the toroidal magnetic field for each 
configuration at high z. The maximum value is close to the disk 
axis in the first two cases, as expected for a collimated wind. The 
last two cases exhibit a much weaker field strength, explaining 
why these flows show no collimation. 



Figure 2. Radial profiles of physical quantities at t = 400. The cuts shown in th upper left panel, (a), were taken at z ~ 0.0 close to the 
disk surface while the cuts in panels (b), (c) and (d) were taken further out at z ~ 72.0. 
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Figure 3. The Alfven lever arm for the four configurations, at 
t = 400, as compared to the prediction from steady-state theory 
(squares). 

The initial magnetic configurations in our studies are 
chosen so that no Lorentz force is exerted on the initial (non- 
rotating) hydrostatic corona described above; more specifi- 
cally, we used initial current-free configurations J = in the 
computational domain. The initial vector potential is sub- 
ject to the boundary condition on the disk that the poloidal 
field on the disk surface (r o ,0) has a power- law form and 
field structure given by equations (23) and (24) respectively. 

An analytic solution for the field throughout the com- 
putational volume can be found for the choice n — on the 
disk. This is the so-called "potential" configuration that we 
used in several of our papers. We can accommodate more 
general initial magnetic configurations however, by using a 
Hankel transform technique to solve for the initial vector 
potential that is subject to these constraints (see Appen- 
dices A and B for details). We chose our particular Hankel 
transform method as the one that best maintained the ini- 
tial equilibrium state for arbitrary amounts of machine time 
(see Appendix B). We always use an unsoftened gravita- 
tional potential from the central object. The initial corona 
is ultimately swept away by the disk wind in all of these 
simulations so that the mass source for our jets becomes 
dominated by the material injected from the fixed disk be- 
low. 

We simulated a representative set of models by choosing 
an initial BP82 model (jj, = -1/4), PP92 (p = -1/2), and a 
yet more steeply declining magnetic field such as /j, — —3/4 
(referred to as the 7Q configuration). These initial config- 
urations are plotted out in the left panels of fig. Q From 
the current profile given in equation (24), we would predict 
that if these flows settle into a stationary state, then the 
models with fi = 0; — 1/4 should be cylindrically collimated 
while the remaining two models; fi = —1/2; —3/4 should be 
"wide-angle" or parabolically collimated flows. This was our 



rationale in choosing these special, illustrative models. The 
mass loading functions corresponding to these four magnetic 
field distributions are progressively less steeply raked; going 
from: 

fc(r )ocrv\ r~ 3/4 , r~ 1/2 , r~ 1/4 (33) 
respectively. 

Our simulations of thin disks require that we specify 
five physical quantities at all points of the disk surface at 
all times. There are the disk density p(r ); components of 
the vertical and toroidal magnetic field, B z (r ) and B$(r )\ 
and velocity components in the disk, v z (r ) and v^,(r ). 
The remaining field component B r (r ) is determined by the 
solenoidal condition while the radial inflow speed through 
the disk v r (r ) ~ for the purposes of the simulations since 
it is far smaller than the sound speed in a real disk. The 
model is described by five parameters, whose values describe 
conditions at the inner edge of the disk at radius r — n. 
Three of these parameters describe the initial corona: /?j 
which is typically 1/3 (and which therefore falls with radius); 
the density jump across the corona/disk boundary r\i which 
is typically 100; and ratio of the Keplerian to thermal energy 
density Si, set to be 300 initially. Two additional parameters 
describe the disk physics. The parameter, Ui, which scales 
the toroidal field in the disk = Vi/r , was found to be 
unimportant. We emphasize that this last initial condition 
plays no role in the subsequent evolution of the simulation - 
even if the initial toroidal field on the disk surface vanishes, 
the wind itself quickly generates the self-consistent field. We 
have never observed that this parameter is important. 

We ran high resolution simulations in which (500 x 200) 
spatial zones were used to resolve a physical region of (80ri x 
20ri) in the z and r directions, respectively. Our simulations 
ran up to 400ti (where U is the Kepler time for an orbit 
at the inner edge of the disk, n). The standard parameter 
settings were 

{•ni,Vi,Vi nj ,5 u [}i) = (100.0,1.0,0.001,300.0,1/3) (34) 

We focused our simulations on the four initial magnetic 
field configurations discussed above - which between them 
give a very broad range of mass load profiles and behaviours 
for the jets. We found that we could only run the BP82 
configuration to t — 400, and barely that, as the simulation 
ground down terribly. Once the density in the corona had 
been cleared out, the Alfven velocity went very high, making 
the timestep prohibitively small. This was also a problem 
with the PP92 case. Since the 7Q case never really blew out 
a lot of material, we managed to run it out to t = 500. 



5 RESULTS 

5.1 Jet collimation - magnetic field lines, density, 
and velocity behaviour 

Figure plots the magnetic field lines of both the initial 
states of each model (left panels) as well as the final states 
of simulations. In addition to the potential case (fi = 0.0), 
simulations of the Blandford-Payne (fj, — —0.25), Pelletier- 
Pudritz (fi = —0.5), as well as an even more steeply raked 
(/x = —0.75) initial magnetic configuration are shown. It is 
evident that the field lines in these right hand panels are 
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(/J = -0.55) 




20 +0 60 



(a) Density structure of jets from the potential (fi = 0, left pan- 
els), and Blandford-Payne (fi = —0.25, right panels) solutions. 
Note the densest material moves along the jet axis. 



progressively less well collimated towards the axis, sugges- 
tive of a wide angle flow at larger radii. A careful examina- 
tion shows that while the potential and BP82 cases achieve 
cylindrical collimation within our computational grid, the 
PP92 and 7Q cases remain more "open". While collimation 
may still occur on spatial scales much larger than our com- 
putational volume, this is difficult to investigate because of 
the very small time step required in regions of large Alfven 
speeds. We see a difference in collimation behaviour that 
agrees with the predictions of steady-state jet theory that 
we outlined in §3.2. It appears that the PP92 case is the 
transition between solutions that collimate to cylinders, and 
those that have parabolic structure in the asymptotic limit. 

Figure 1 also shows the positions of the Alfven surfaces 
for each configuration marked on the field lines. The BP82 
final state shows that the Alfven surface is indeed nearly 
a cone as expected while the PP92 simulation qualitatively 
agrees with the expected scaling ta tx rj/ 2 (see §3.2). The 
details of the outflow for the "potential" configuration are 
found in OPS97 and OP97a. 

The mass load function for these models is shown in 
Figure |2(a)| wherein the values of the function k(r,z) are 
numerically tabulated along a radial spatial cut through the 
flow that is a few pixels above the disk inside the computa- 
tional zone. In stationary state, the mass loads should settle 
down to power-law radial distributions that match the input 
at the surface of the disk, predicted in equation (31). The 
figure shows power-law fits for physical quantities tj) (in this 
case, ip = k) of the form 



Ijj. = -0.5) 
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(b) The corona has been cleared out by the uncollimatcd disk 
wind. The 'jets' from the Pelletier-Pudritz (/i = —0.5, left pan- 
els), and steeper (/i = —0.75, right panels) solutions should be 
compared to the velocity profiles in Fig. |3] 

Figure 4. Contours mark surfaces of constant density in the jet 
at four different times (t= 60, 120, 200, 400). 



which give for the 4 different magnetic configurations (al- 
ways starting with the potential configuration fj, = 0), 
a = —0.98, —0.78, —0.68, —0.57, which compares favourably 
with the predictions except for the last case, the 7Q model. 
The simulations of the potential flow and the BP82 outflows 
have settled into stationary states so that the mass loads in 
the jet closely give the input value at the disk surface (-1.0 
and -0.75 respectively). The PP92 simulation has not set- 
tled down completely by the end of our run, but well enough 
that there is a good correspondence with the predicted load 
profile of -0.5 compared to our numerical value -0.68. The 
greatest discrepancy is found for the 7Q run, but as noted, 
this has not achieved a completely stationary state even by 
the end of t = 400. 

A quantitative analysis of the radial density behaviour 
of these 4 outflows is given in Figure [2(b) | The potential and 
BP82 configurations have greater peak densities than the 
PP92 and 7Q models. Another difference appears to be that 
whereas the former two have a declining power-law density 
behaviour with 



-0.74 (p = 0);a = -0.34 {(j, = -1/4) 



(36) 



- the potential and BP82 cases respectively - the radial be- 
haviour is quite different in the PP92 and 7Q cases which 
appear to be nearly independent of radius. In every case, 
there is a central "spike" in the gas density which is the 
material at the smallest radius that moves more parallel to 
the outflow axis. The radial position of this "spike" moves 
outward in radius as one can see from the figures. 

The corresponding radial behaviour of the poloidal ve- 
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locity is analyzed in Figure [2(c) | by means of a spatial cut 
through the flow taken most of the way down the jet at 
z = 72.08. The peak speed in each of the first 3 models in- 
creases. There is also a distinct fall-off in jet outflow speed 
with disk radius, which for the potential and BP82 cases 
behave as 

a = -0.78, a = -0.81; (37) 

which are virtually identical. It is evident from this figure 
that the PP92 and 7Q cases have a larger "hole" in the 
middle of their jets, and that their poloidal speeds fall off 
less steeply with jet radius. 

The radial profile of the toroidal magnetic field for the 
models is shown in Figure [2 (d)| In all of our simulations, the 
toroidal field must vanish on the rotation axis. It reaches a 
peak value at the inner edge of the outflow and then de- 
clines. This means that the toroidal magnetic field pressure 
exerts an inward pressure force from this peak, and an out- 
wards pressure gradient force at larger radii than the peak 
radius (see OP97a). While the potential and BP82 cases 
have similar radial power law behaviour for the toroidal field; 
B^ir) oc r -0 ' 64 on average at the larger radii, the PP92 and 
7Q cases have much weaker toroidal fields everywhere. 

We compare the predicted structure of the Alfven sur- 
face and that from theory in Figure 3 where we plot the 
Alfven lever arm versus the footpoint radius of a given field 
line. The analysis of equation (25) predicts that as one goes 
from jit = — ► = —3/4, then 

r A /r <xrl /2 , const., r~ 1/2 , r' 1 (38) 

respectively. The first three solutions compare well to the 
theory. The steepest case has its Alfven surface too close 
to the disk surface to be well-resolved, which may explain 
the discrepancy. We see that the Alfven lever arm usually 
lies within the domain of 2-3 in value for the innermost field 
lines of the flow. 

The evolution of the density structure of these jets is 
shown in the snapshots presented in Figures |4(a)[|4(b)| for 
4 different times (t = 60, 120,200,400). In each simulation, 
we see the jet driving a bow shock through the quiescent 
gas in the first frame. In the second frame the shock has 
just moved out of the grid. In all four of these simulations, 
one clearly sees that the densest gas in the jet is always well 
collimated and moving nearly parallel to the axis. This is 
a general consequence of any reasonable mass loading for a 
jet. Thus, whether or not the outflow has a wide-angle as- 
pect to it or not, the bulk of the material that we see in 
all of our simulations moves along in a fairly concentrated 
and collimated manner. Thus, all of our simulations when 
observed in appropriate forbidde n line emission, w ould ap- 
pear "jet-like" - as first noted bv lShu et all §000) for their 
X-wind simulations. 

The poloidal velocity vectors associated with these out- 
flows is shown in Figures |5(a)H5(d")| The four frames shown 
in each of these panels correspond to the 4 times shown in 
Figs. |4(a)|4(b)| A well-collimated jet is visible in the velocity 
figures for the first two cases - potential and BP82 outflows. 
In the fourth case (7Q) we see that whereas the jet density 
is highly collimated (Fig. 4b, right panel) along the axis - a 
wide-angle velocity profile is obvious. Thus, the bulk of the 
jet material follows the few field lines that have collimated 
parallel to the disk axis, but is clearly sub-Alfvenic. The 



bulk of the kinetic energy of the flow is in the wide-angle 
wind. Thus, the PP92 and 7Q cases represent two models in 
which a wide-angle component could carry momentum suf- 
ficient to produce a wide-angle molecular outflow. The final 
stage (t = 400) of the PP92 solution resembles the 7Q case 
more closely than either of the well-collimated solutions. 

Taken together, Figures show the existence of two 
regimes of disk-driven winds; one that results in a highly col- 
limated jet, and the other that produces a wide-angle wind. 
The PP92 configuration seems to be a transitional model 
lying between these two distinct regimes. In all situations, 
the bulk of the flow's mass density is fairly well collimated 
along the flow axis. The cylindrical or parabolic collimations 
(/J. = 0,-1/4 and /i = —1/2,-3/4 cases respectively) pri- 
marily affect the velocity field and a small amount of the 
gas in the outflow. 

5.2 Jet dynamics - moving along a field line 

In Figs. |6(a)| — |6(d)| we trace the behaviour of physical 
quantities along the fifth innermost fieldline which is chosen 
since it remains on the computational mesh and is concen- 
trated more towards the central axis of the flow. Comparing 
the upper-left panel for each configuration, we see that the 
potential case reaches an Alfven Mach number of about 1.5, 
but never reaches its fast magnetosonic speed. The BP82 
case accelerates much more quickly, and has reached its fast 
magnetosonic point at z — 1.5, reaching a peak Alfven Mach 
number of 7 before the poloidal velocity levels off. 

The PP92 configuration reaches a very high Alfven 
Mach number of nearly 50 before it undergoes an abrupt 
downturn in value. The 7Q case does not really seem to 
settle into as stationary a configuration, but the results do 
seem to suggest that it is in a different regime of Alfven 
Mach number. 

The rather abrupt decrease in the outflow Mach num- 
bers seen in this figure arises from the presence of shocks and 
knots along the innermost field lines. This is evident from an 
examination of the Alfven and FM points along inner field 
lines shown in Figure Q While the potential flow seems to 
have reached a very stable, stationary state by the end of 
the simulation, the BP82 and PP92 simulations require in- 
creasingly longer times to settle into stationary flows while 
the 7Q case did not settle into a stationary state in the far 
field region at z/r% 60. Another indication that shocks are 
responsible for these abrupt changes is that the toroidal to 
poloidal magnetic field, as well as the value of the toroidal 
flow velocity, also undergo similar sharp changes in value. 

The outflow velocity along the field line, in units of the 
Kepler speed at the foot of the field line, increases for these 
first three models. Thus the highest speed flow in the centre- 
most region of a jet arises in the PP92 model, more even than 
the BP82 solution. As for the rotation speed of material in 
the outflow, it always falls to a fraction of its initial value, as 
theory predicts (because of the expansion of the flow, as well 
as the transport of angular momentum that is carried by the 
twisted field) . We defer our analysis to the radial rotational 
profile of material in the jet to §5.4 below. 
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20 40 GO SO 

'A, 

(a) Potential case: Note that the highest velocities are nearest the outflow axis, and that the flow 
is very well-collimated 



fJ, = — 0.25 time = 60.0 




20 40 GO 80 



(b) BP82 case: Note that the highest velocities are nearest the outflow axis, and that the flow 
gently recollimates 

Figure 5. Poloidal velocity fields of the 4 models, shown at the 4 different times in Fig. |2(a)] 
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/j. = -0.5 




time - 200.0 



time = 400.0 



(c) PP92 case: The flow seems to be collimating in the early stages of the jet evolution, but as can 
be seen in the last panel, the flow assumes a conical geometry. 



H = -0.75 



time = 60.0 



time - 120.0 



time - 200.0 




(d) 7Q case: The steepest configuration barely collimates at all and evolves differently from the 
other solutions. The wide-angle wind is even more extreme than the Pelletier-Pudritz solution 



Figure 5. ...continued. 
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Figure 6. Physical quantities along the fifth innermost fieldline, plotted against z/r. The position of the Alfven and fast magnctosonic 
(marked by FM) critical points are shown. 
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Figure 7. Quantification of the energies involved in the flow. In the upper/lower panel, the poloidal magnetic/kinetic energy is compared 
to the toroidal part. Note that the the poloidal magnetic energy in the last case is double that in the other cases. 
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Finally, the ratio of the toroidal to poloidal magnetic 
field along our illustrative field line clearly shows that the 
predominant magnetic field in jets beyond their FM surfaces 
is the toroidal field component. The stability of a jet that is 
dominated by a toroidal field is an important problem, and 
the 3D simulations needed to investigate the importance of 
kink modes are still rare. Our own 3D study of a jet in 
an initially uniform field JOuved. Clarke, fc PudriO [2003 ) 
shows that jets survive the threatening non-axisymmetric 
kink (m = 1) mode, and we plan to extend this work to 
these more general configurations. 



5.3 Jet energetics 

The bulk energetics of our 4 model outflows, as a function of 
time, are shown in |7(a)|j7(d)| In every case, we see that the 
jets end up being dominated by the poloidal kinetic energy of 
the jet, followed in magnitude by the energy in the poloidal 
magnetic field, and then the energy in the toroidal magnetic 
field. The lowest energy in all of the outflows is the energy 
in the overall rotation of the jet. Thus, we find that 



E, 



kin,pol 



mag,pol 



> E„ 



> E k 



(39) 



These energies shown in this figure are given in units of 
the following physical quantities: the poloidal and toroidal 
kinetic energies are in units of 2-Kp i v\ i rj, while the 
poloidal and toroidal magnetic energies are given in units 
of (BlJ2)rl 

The ratios of the various energy terms are more interest- 
ing. In going through the 4 cases /i — to fi = —3/4, these 
figures show that the ratio of the poloidal kinetic energy 
(bulk energy in outflow) to the bulk energy in the poloidal 
fields, takes the values: 



5.4 Jet rotation 

The rotational velocity (left hand column) and toroidal mag- 
netic field strength (right hand column) throughout the 4 
outflows, are shown in Figure |H| As was found in the po- 
tential configurations of OPS97, the top panel in this figure 
shows a beautifully, well collimated outflow with the toroidal 
field components nearly constant along cylinders. This is 
also seen for the BP82 case although, not as strongly as the 
potential case. The wide-angle behaviour of the final two 
cases is apparent in both the toroidal velocity and field con- 
tours as well. The explanation for this wide-angle behaviour 
lies in the distribution and strength of the toroidal magnetic 
field (right panels), which is too weak to collimate the flow, 
as contrasted with the upper two cases, which clearly has a 
strong £>0 along the jet axis. 

The radial profile of the rotation speed of material in 
the outflows are shown in Figure [0] The two columns show 
cuts at different distances z down the jet axis, the right hand 
column being within the Alfven surface, while the left hand 
column is at high z. At high z, the outermost radial profiles 
of the outflow's rotation speed for the first 3 models can be 
well fit by power laws; 



v<f,(r) oc r a 



-0.76 



0.66; 



-0.46 



(42) 



(43) 



respectively, whereas the 7Q simulation is even flatter in 
profile but too variable to fit with a powerlaw at t = 400. 

These rotational profiles are sufficiently different that it 
may be possible for future observations to place constraints 
on the precise nature of the mass load by mapping out the 
radial rotational profile, and using it to infer the underlying 
mass load. 



Eki 



E 



mag,pol 



= 1.5, 2.0, 



6.3. 



(40) 



This shows that the most efficient acceleration of the flows 
and the conversion of gravitational binding energy into bulk 
flow is from the less steeply raked mass loads - ie the flows 
that carry the most mass. We noted earlier that the PP92 
solutions achieved the highest terminal Alfven Mach num- 
bers, which goes together with this result. 

A comparison of the energy in the poloidal and toroidal 
field components is most revealing: 



f-^mag ,pol 
Emag ,tor 



= 1.2, 1.3, 2.5, 8.0. 



(41) 



Clearly, the relative importance of the toroidal field be- 
comes progressively less as one goes through this mass load 
sequence - exactly as our theoretical model predicts. The 
potential and BP82 solutions (fi — 0, —0.25) both have 
toroidal magnetic energies comparable to their poloidal en- 
ergies, whereas the other two solutions are dominated by 
the poloidal field component. While all four configurations 
have comparable toroidal magnetic energies (with the colli- 
mated solutions having more at later times), the energy of 
the poloidal magnetic field in the last case is almost three 
times larger than the others. 



6 DISCUSSION AND CONCLUSIONS 

We have demonstrated that accretion disks play a major 
role in controlling the collimation of disk winds through the 
radial mass load profiles that they impose on these outflows. 
The conservation laws for stationary, axisymmetric outflows 
show that the toroidal field at any point in the jet is reg- 
ulated by this mass load function. Given that the toroidal 
field is responsible for the collimating hoop stress in the jet, 
accretion disks can help determine the structure of jets and 
outflows far from their point of origin on the disk. 

Our simulations divide into two classes - those that 
achieve collimation towards a cylinder, and those that have 
a wide-angle structure. We showed that these correspond 
quite well to predictions made by HN89 about jet collima- 
tion, namely, that if the current I(r) goes to zero in the 
limit of large jet radius, that such field lines would have a 
parabolic structure, whereas currents that remain finite or 
diverge result in jets that collimate cylindrically. 

We argue that the differences between wide-angle and 
highly collimated jets and outflows inferred in the obser- 
vational literature, are two aspects of the same underly- 
ing mechanism. Hydromagnetic disk winds can achieve both 
types of configuration depending on the mass loading that 
takes place in the underlying accretion disk. Specifically, 
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Figure 8. Contours of the toroidal velocity (left panels) and toroidal magnetic field (right panels) at t = 400 for the four cases. This 
figure can be compared to |5(a)|5(d)| Note that the magnetic field in the last two panels is too weak to collimate the flow. 



we find that the potential and Blandford-Payne magnetic 
field distributions B p (r ) oc r" 1 , To'^ A respectively, result in 
mass loading profiles that produce cylindrically collimated 
jets. The Pelletier-Pudritz configuration on the other hand, 
B p (r ) oc r ( 7 3/ ' 2 results in a finite current and is the bound- 
ary between cylindrically and parabolically collimated out- 
flows. It may be that nature produces magnetic configu- 
rations in accretion disks that arise through an optimiza- 
tion process (eg. minimum magnetic energy configuration) 
of some kind. We speculate that it may be that small fluctu- 
ations in jet mass loading, about an energetically preferred 
PP92 profile, would result in jets that are either cylindri- 
cally, or parabollically (ie wide-angle) collimated. This could 
explain the variety of jet collimations that are suggested by 
the observations. 

Our res ults allow one to understand why simulation of 
X- winds by IShu et ail |2000) , or split monopoles by Ro- 
manova et al (1997) result in wide-angle flows whose field 
lines do not collimate towards the outflow axis. This is a 
consequence of having steep magnetic field gradients, and 
hence shallow mass loading profiles which, as we have seen, 
result in weak toroidal fields in the jet. T hese are unable to 
collimate the nearly radial field lines. As IShu et alJ (l200fj) 
note, this does not prevent the densest material in the out- 



flow from concentrating in flow that is parallel to the outflow 
axis. We see this in all of our simulations. 

The radial profiles of various outflow quantities such as 
the poloidal outflow speed v po i are tabulated in this paper. 
Our outflows all show the "onion- layer" velocity structure 
that are revealed by the observations of protostellar jets, 
namely, that the highest velocity gas is in the interior of the 
jet while the lower velocity material are found at increasingly 
larger outflow radii. Of particular interest is the behaviour 
of jet rotation v$(r) oc r a as one goes radially through the 
jet. We find that this distribution depends on the underlying 
mass load (and hence magnetic configuration in these simu- 
lations). Specifically, a — —0.76, —0.66, —0.46 as one goes 
from the potential configuration, to the Blandford-Payne 
distribution, and on to the Pelletier-Pudritz distribution of 
disk fields /i = 0, —1/4, —1/2 respectively. 

In conclusion, this paper shows that accretion disks can 
exert a long-range control of jet properties such as their 
collimation and spin, as a consequence of their mass loading. 
What remains to be clarified is exactly how the accretion 
disk actually does this. It may be that the accretion disk 
corona plays an important role in this process. Dynamical 
3D global simulations of the disk and its outflow, as now 
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Figure 9. Radial profiles of the toroidal velocity within the Alfven surface (left panels) and at high z (right panels) for each solution at 
t = 400. The sub-Alfvenic flow can easily be fit by a power law. At high z, only the outermost radii can be well fit, while the steepest 
case still exhibits complex behaviour. 



being carried out by an increasing number of authors, may 
offer the best way of exploring the physics of this problem. 
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APPENDIX A: COMPUTING THE INITIAL 
CONFIGURATION OF THE MAGNETIC 
FIELDS 

First, we assume axisymmetry, and thus all -J=^ terms arc 
ignored. Second, we are considering only poloidal compo- 
nents; all toroidal magnetic field arises from the centrifugal 
motion of the disk. If we take the magnetic field as 



B = V x A 



then 



and 



' dz ' 

=0; 
ldrA ± 



(Al) 



(A2) 



(A3) 



(A4) 



r dr 

Thus it remains to find some function A^ which satis- 
fies our boundary and initial conditions. We considered only 
solutions that were current-free, thus keeping the Lorentz 
force F = J x B on the object also equal to zero. Since 
J = V x B, this means the current-free condition can be 
written 

dB r _ dB z 
dz dr 
which gives us a separable equation for 

& . , 

(A6) 



d 2 A 4> I d A* _ Aj, 
dr 2 r dr r 2 







9 A * = -k 2 



(A5) 



dz 2 



with (non-divergent) basis functions A^^ 
A^fc = Ji(fer) exp(— k\z\) 



(A7) 
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yielding a general class of solutions 

P oo 

A 4> (r,z) = S(fc)Ji(A;r)exp(-fc|z|)dfe (A8) 
Jo 

It now remains to find the amplitude S(k), from some 
boundary condition. We exploit this freedom and let B z in 
the disk have a radial dependence given by: 

B z (r,z = 0) = br"- 1 (A9) 

where b is a normalizing constant. 

Using Eqn A. 4 and A. 8 and exploiting the properties of 
Bessel functions, we have the following form for B z : 

P oo 

B z (r,z = 0)= S(k)J (kr)kdk (A10) 
Jo 

The inverse Hankel Transform gives us an equation for S(k) 

p oc 

S(k) = / B z {r,z = 0)J o (kr)rdr 

(All) 



= b r M Jo(kr)dr 
Jo 



which has the solution 



where 



r(i/2 + M /2) 



/(A*) = 2 



(A12) 



(A13) 



r(l/2- M /2) 

(Note that this analytic transform was useful for testing our 
numerical integrator, as explained in the next section.) 
Thus, equation A. 8 becomes 

POO , 

Mr,z)=bf(ji)J ^ rTT J 1 (fcr)exp(-fc|^|)dfc (A14) 

This Hankel transform has only two analytical solu- 
tions, when (i — +1, which has vertical field lines (B z — b, 
where b is a constant), and when n — 0, which is the poten- 
tial solution discussed in OPI. 

We employed a numeric integrator to evaluate three 
additional configurations, approaching the lower limit of 
/i = — 1; jU = —1/4, which corresponds to the self-similar so- 
lutions of BP82, /i = —1/2, the solution suggested by PP92, 
and fj, = —3/4. 



and 



well near the analytic limits (// = — 1, +1) of the functions we 
were evaluating. It was capable of integrating Bessel func- 
tions of any order, but since our transforms were only of 
order and 1, thi s feature was unn ecessary. Anderson's pro- 
gram 'ZHANKS' dAndersonlll97sl) is a single precision code 
that uses adaptive filter weights. The filters are designed by 
using Hankel Transforms with known analytic solutions: 



Aexp(-aA 2 )J (bA)dA = [exp(-6 2 /4a)]/(2a) (Bl) 



A exp(— aA ) Ji (b\)d\ = &[exp(-b 2 /4a)]/(2a) 2 (B2) 



where a > 0, b > 0. 

Its algorithm uses the fact that algebraically related ker- 
nel functions (the function being transformed by convolution 
with the Bessel function) transform the same way. Thus, the 
computed kernal function is saved on its first call, and subse- 
quent calculations use the results of previous calculations to 
evaluate the transforms, making computation much faster. 
The solutions computed by the ZHANKS algorithm are also 
behaved better at the analytic limits, matching closely with 
the solutions we tested using MATLAB. 

In addition to comparisons with known analytic solu- 
tions, the results of the algorithms were tested in a hydro- 
static run of our simulation in which the disk did not rotate. 
Since the corona was in pressure balance with the central 
star and the surrounding disk, only magnetic forces would 
disturb the gas. As the initial conditions should be com- 
pletely current-free, no Lorentz forces should be developed, 
and the gas should remain motionless. This was in fact the 
case, except for slight motions arising from a combination 
of the accuracy of the algorithms and the discrete nature of 
the grid. Our hydrostatic simulations were run out beyond 
400 timesteps for each magnetic configuration. The motions 
were slightly increased for the more radially dependent mag- 
netic profiles (fi = —1/2 and fi — —3/4)), but in all cases 
their movement was such that the hydrostatic balance was 
disturbed to no more than 5ri after 400£,: 



APPENDIX B: HANKEL TRANSFORM 
ALGORITHM 

Since there are a bevy of Hankel transform algorithms avail- 
able in the public domain, it was unnecessary to devise our 
own. In our search for a numeric integration program ca- 
pable of performing Hankel transforms, we tested several 
ap plications. The two most te s ted p rograms were designed 
bv lAndersonl dl979T l: IPiessens! (I1982T) . The most important 
criteria for our choice were speed and accuracy. The initial 
conditions for our 500x200 grid were set point by point, so 
it was important that the transforms be done quickly and 
efficiently. 

Piessens' program 'hankel' is a double-precision code 
using Fast Fourier Transforms (which are done by gaus- 
sian quadrature methods) to evaluate each transform. It ran 
somewhat slowly, and the solutions it gave did not behave 



